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Abstract 

Abstract.  We  prove  mathematically  that  in  order  to  avoid  point-optimization  at  the  sampled  design  points 
for  multipoint  airfoil  optimization,  the  number  of  design  points  must  be  greater  than  the  number  of 
free-design  variables.  To  overcome  point-optimization  at  the  sampled  design  points,  a  robust  airfoil 
optimiza-  tion  method  (called  the  pro  le  optimization  method)  is  developed  and  analyzed.  This 
optimization  method  aims  at  a  consistent  drag  reduction  over  a  given  Mach  range  and  has  three 
advantages:  (a)  it  prevents  severe  degradation  in  the  o-design  performance  by  using  a  smart  descent 
direction  in  each  optimization  iteration,  (b)  there  is  no  random  airfoil  shape  distortion  for  any  iterate  it 
generates,  and  (c)  it  allows  a  designer  to  make  a  trade-o  between  a  truly  optimized  airfoil  and  the  amount 
of  computing  time  consumed.  For  illustration  purposes,  we  use  the  pro  le  optimization  method  to  solve  a 
lift-constrained  drag  minimization  problem  for  2-D  airfoil  in  Euler  ow  with  20  free-design  variables.  A 
comparison  with  other  airfoil  optimization  methods  is  also  included. 
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Abstract.  We  prove  mathematically  that  in  order  to  avoid  point-optimization  at  the  sampled  design 
points  for  multipoint  airfoil  optimization,  the  number  of  design  points  must  be  greater  than  the  number  of 
free-design  variables.  To  overcome  point-optimization  at  the  sampled  design  points,  a  robust  airfoil  optimiza¬ 
tion  method  (called  the  profile  optimization  method)  is  developed  and  analyzed.  This  optimization  method 
aims  at  a  consistent  drag  reduction  over  a  given  Mach  range  and  has  three  advantages:  (a)  it  prevents  severe 
degradation  in  the  off- design  performance  by  using  a  smart  descent  direction  in  each  optimization  iteration, 
(b)  there  is  no  random  airfoil  shape  distortion  for  any  iterate  it  generates,  and  (c)  it  allows  a  designer  to  make 
a  trade-off  between  a  truly  optimized  airfoil  and  the  amount  of  computing  time  consumed.  For  illustration 
purposes,  we  use  the  profile  optimization  method  to  solve  a  lift-constrained  drag  minimization  problem  for 
2-D  airfoil  in  Euler  flow  with  20  free-design  variables.  A  comparison  with  other  airfoil  optimization  methods 
is  also  included. 
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Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  Optimization  of  aerospace  vehicles  or  aircraft  wings  is  based  on  a  mathematical 
model  of  the  physical  reality.  The  design  variables  are  parameters  in  the  mathematical  model  and  changes 
in  these  design  variables  result  in  new  physical  objects.  Aerodynamic  optimization  is  a  process  of  finding 
a  set  of  design  variables  that  corresponds  to  a  new  physical  object  with  better  aerodynamic  and  structural 
properties. 

For  airfoil  shape  optimization,  design  variables  are  parameters  that  define  the  airfoil  shape.  One  accepted 
practice  for  modeling  airfoil  shapes  is  to  use  empirical  algebraic  expressions  that  are  based  on  the  knowledge 
of  aerodynamic  properties  of  airfoils.  There  are  two  advantages  to  such  an  airfoil  model: 

(a)  Each  set  of  design  variables  (such  as  maximum  thickness,  camber,  radius  of  leading  edge,  etc.) 
generates  an  airfoil  with  aerodynamic  properties  that  are  well  understood  by  experts. 

(b)  The  number  of  design  variables  is  small,  thus  the  corresponding  optimization  problem  is  computa¬ 
tionally  affordable. 

A  drawback  to  such  an  airfoil  model  is  that  the  optimization  process  is  merely  selecting  a  desirable  airfoil 
shape  among  the  predetermined  shapes  given  by  the  specific  airfoil  model.  The  usefulness  of  the  optimal 
airfoil  is  essentially  determined  by  the  usefulness  of  the  airfoil  model. 

To  achieve  truly  innovative  airfoil  designs,  it  is  therefore  desirable  to  consider  “free^forrn”  shape  opti¬ 
mization,  which  allows  “truly  optimum”  designs  to  be  computed  [5,  Section  1].  Here  “free-forrn”  means 
geometric  shapes  represented  by  a  linear  combination  of  general  basis  functions  (such  as  splines  or  sinusoidal 
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basis  functions).  However,  there  are  several  challenges  associated  with  ‘‘free-forrn”  airfoil  shape  optimiza¬ 
tion.  Drela  [5,  Subsection  5.1]  pointed  out  that  if  presented  with  sufficient  design  model  resolution,  an 
optimizer  will  readily  (and  annoyingly)  manipulate  and  exploit  the  flow  at  the  smallest  significant  physical 
scales  present  that  tends  to  produce  improved  performance  only  near  the  sampled  operating  conditions. 
The  point-optimized  airfoil  often  shows  a  possibly  severe  degradation  in  the  off-design  performance  and 
optimized  aerodynamic  shapes  are  usually  ‘hioisy”  and  usually  require  a  posteriori  smoothing. 

The  main  objective  of  this  paper  is  to  develop  a  robust  airfoil  optimization  scheme  for  achieving  consistent 
drag  reduction  over  a  given  Mach  range.  This  scheme  (called  the  profile  optimization  method)  has  the 
following  three  advantages:  (a)  it  prevents  severe  degradation  in  the  off-design  performance  by  using  a  smart 
descent  direction  in  each  optimization  iteration,  (b)  there  is  no  random  airfoil  shape  distortion  for  any  iterate 
it  generates,  and  (c)  it  allows  a  designer  to  make  a  trade-off  between  a  truly  optimized  airfoil  and  the  amount 
of  computing  time  consumed. 

The  term  ‘Robustness”  has  been  used  to  mean  a  variety  of  things.  For  optimization  under  uncertainty, 
we  can  distinguish  the  following  meanings  and  goals  of  “robust  optimization” : 

1.  Identify  designs  that  minimize  the  variability  of  the  performance  under  uncertain  operating  condi¬ 
tions.  This  is  the  objective  of  Taguchi  methods  [8],  which  are  most  practical  when  the  expected 
value  of  the  performance  can  be  adjusted  at  negligible  cost. 

2.  Mitigate  the  detrimental  effects  of  the  worst-case  performance.  This  is  the  objective  of  rninirnax 
strategies,  which  can  be  used  to  find  a  design  with  the  optimal  worst-case  performance  [4]. 

3.  Provide  the  best  overall  performance  of  a  system  by  maximizing  the  expected  value  of  its  utility 

[12]. 

4.  Achieve  consistent  improvements  of  the  performance  over  a  given  range  of  uncertainty  parameters. 
This  is  the  main  objective  of  this  work. 

The  paper  is  organized  as  follows.  In  Section  2,  we  give  a  brief  introduction  to  the  airfoil  optimization 
problem.  Section  3  is  devoted  to  Prelaws  hypothesis  on  the  necessity  of  using  at  least  0{m)  design  points 
for  multipoint  airfoil  shape  optimization,  where  rri  is  the  number  of  free-design  variables.  The  main  result 
is  a  mathematical  proof  of  the  fact  that  the  multipoint  airfoil  shape  optimization  needs  at  least  (m  +  1) 
design  points.  We  present  two  robust  airfoil  shape  optimization  formulations  in  Section  4  and  prove  the 
mathematical  equivalence  between  these  two  formulations  under  finite  sampling  in  Section  5.  The  profile 
optimization  method  is  introduced  in  Section  6.  The  implementation  issues  of  the  profile  optimization 
method  are  discussed  in  Section  7  and  numerical  simulation  results  are  given  in  Section  8.  Final  conclusions 
are  drawn  in  Section  9. 

2.  Airfoil  Shape  Optimization.  Recently,  there  has  been  significant  progress  in  airfoil  shape  opti¬ 
mization  (see  [2,  3,  5,  16]  and  references  therein).  These  papers  demonstrate  impressive  shape  optimization 
using  high-fidelity  CFD  codes,  reliable  grid  generation,  and  numerically  efficient  sensitivity  calculations. 
Equally  impressive  progress  has  been  made  in  optimization  of  3-D  wings  [6,  7,  17]  and  in  coupled  structural- 
aerodynamic  optimization  [9].  However,  except  Drela’s  work  [5],  these  aerodynamic  shape  optimization 
projects  all  find  optimal  shapes  based  on  fixed  operating  conditions.  In  this  paper,  we  study  a  simplified 
2-D  airfoil  shape  optimization  problem  using  a  low-fidelity  Euler  flow  solver,  but  we  include  uncertainty  in 
the  operating  conditions  for  the  airfoil  shape  optimization. 

Airfoil  shape  optimization  is  a  PD  E-const  rained  optimization  problem.  A  general  mathematical  frame¬ 
work  for  airfoil  optimization  can  be  described  as  follows  (see  [3,  16,  17]  for  details). 

For  a  given  flow  and  turbulence  model  (a  set  of  PDEs),  one  could  compute  the  energy,  velocity,  and 
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pressure  around  a  given  airfoil,  which  will  be  denoted  by  one  vector-valued  function  S,  For  an  airfoil 
optimization  problem  (which  is  a  2-D  problem),  5  is  a  function  of  {x%y)  for  fixed  (£>,a,  Af),  where  D 
denotes  the  set  of  geometric  design  parameters  that  defines  the  airfoil,  a  is  the  angle  of  attack,  and  M  is 
a  given  Mach  number.  Note  that  M  and  a  are  parameters  representing  the  flight  conditions.  The  PDE 
equations  for  S  are  solved  by  a  numerical  method  that  involves  grid  generation,  discretization  of  the  PDEs, 
and  iterations  for  solving  systems  of  nonlinear  equations. 

For  a  given  geometric  design  vector  I>,  the  generated  grid  (denoted  by  X)  is  determined  by  I>,  i.e., 
X  =  X{D)  is  a  function  of  D,  Usually,  X{D^)  is  generated  for  an  initial  design  vector  based  on  a  grid 
generation  method,  and  X {D)  (for  D  ^  D^)  is  generated  by  a  grid  movement  strategy  (such  as  the  spring 
analogy  or  the  elasticity  model)  [16]. 

The  system  of  nonlinear  equations  derived  from  a  discretization  of  the  PDEs  for  the  flow  and  turbulence 
model  will  be  denoted  by 

(2.1)  R{D,a,M,X,S)  =0, 

where  5  is  a  vector  whose  components  are  values  of  S  at  the  grid  points  (or  the  cell  centers).  For  any  given 
(D,a,  Af,  X),  a  numerical  solution  of  (2.1)  gives  S, 

If  we  consider  minimization  of  the  drag  coefficient  under  a  minimal  requirement  on  the  lift  coefficient, 
then  the  multipoint  airfoil  optimization  problem  can  be  formulated  as  follows  [5]: 

r 

(2.2)  min  E  Wi  Ca (^D,  ai,  Mi,X {D),Si{D, 
subject  to 

(2.3)  ci{^D^ai^Mi^X{D)^Si{D^ai)^  >  for  1  <  i  <  r 

and  D  ^  Xy  where  Wi^s  are  positive  weights,  is  the  minimal  lift  value,  is  a  given  feasible  set  for 
geometric  design  variables,  Si{D^ai)  is  the  solution  of  (2.1)  for  Af  =  Mi  and  a  =  a^,  and  (or  q)  denotes 
an  approximate  value  of  the  drag  (or  lift)  coefficient.  The  feasible  set  X  mainly  depends  on  the  airfoil  model. 
In  this  paper,  we  use  splines  for  airfoil  shape  modeling:  the  components  of  D  are  the  control  points  for  a 
spline  curve  and  the  feasible  set  X  =  ]R^,  where  m  is  the  number  of  geometric  design  variables. 

Recently,  Drela  [5]  studied  the  behavior  of  the  optimization  solutions  of  (2.2)  in  two-dimensional  viscous 
flow  when  the  number  of  free-design  variables  is  relatively  large.  Drela  [5]  concluded  that  increasing  the 
number  of  geometric  design  variables  requires  a  corresponding  increase  in  the  number  of  design  conditions 
(Mach  numbers)  used  in  the  multipoint  optimization  problem  (2.2).  He  also  suggested  that  it  is  necessary  to 
have  r  =  (9(m),  where  m  is  the  number  of  free-design  variables,  to  achieve  a  smooth  airfoil  geometry.  Other 
notable  conclusions  made  by  Drela  [5,  Conclusion]  are: 

•  Near-continuous  sampling  of  the  operating  space  (i.e.,  in  the  range  of  Mach  numbers)  may  be 
required  in  the  theoretical  limit  of  a  general  airfoil  design  problem  with  a  very  large  number  of 
degree  of  freedom  (for  geometric  variables)  —  a  very  expensive  proposition. 

•  The  most  suitable  operating  points  to  be  actually  sampled  in  multipoint  airfoil  optimization  (i.e., 
Afi,  Af2, . . . ,  Afy.)  are  not  apparent  a  priori.  From  limited  experience,  sampling  somewhat  beyond 
the  expected  operating  range  appears  to  be  best. 

•  The  point  weights  (i.e.,  •  •  •  j^r)  used  in  multipoint  airfoil  optimization  are  arbitrary,  and 

their  appropriate  values  can  not  be  easily  estimated  without  prior  experience. 

•  Optimized  aerodynamic  shapes  are  usually  ‘hioisy”  and  require  a  posteriori  smoothing. 
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3.  The  Critical  Number  of  Design  Points  for  Multipoint  Airfoil  Optimization.  In  this  sec¬ 
tion,  we  give  a  mathematical  proof  of  Drela’s  hypothesis  on  the  necessity  of  having  the  number  of  design 
points  proportional  to  the  number  of  design  variables.  Specifically,  we  prove  that  in  order  to  avoid  point- 
optimization  at  the  design  points,  the  number  of  design  points  must  be  at  least  (m  +  1),  where  rri  is  the 
number  of  free- design  variables. 

In  (2.2),  the  angle  of  attack,  a^,  is  predominantly  determined  by  the  constraint  on  the  lift  coefficient  at 
Mi,  Therefore,  it  is  theoretically  possible  to  eliminate  the  constraint  and  consider  the  following  unconstrained 
reformulation  of  (2.2): 

r 

(3.1)  nun  ^Wif{D,Mi), 

where  Af^’s  are  design  points,  are  positive  weights,  and  /  is  a  function  related  to  the  drag  coefficient. 
Let  D  be  an  optimal  solution  of  (3.1).  Then 


(3.2) 


where  ^  denotes  the  gradient  of  /  with  respect  to  D, 

If  r  <  m,  then  (r  —  1)  <  m.  So  we  can  find  a  nonzero  vector  AD  in  the  orthogonal  complement  of  the 
subspace  of  generated  by  the  following  (r  —  1)  vectors: 

. 

By  (3.2)  and  wi  >  0,  AD  must  also  be  orthogonal  to  ^{D^  Mi),  Therefore,  we  obtain  a  nonzero  vector 
AD  such  that 


(3.3) 


Mi)^  AD^  =0  for  1  <  i  <  r. 


where  {u^v)  denotes  the  dot  product  of  vectors  u  and  v. 
Let  be  any  Mach  number  such  that 


(3.4) 


Let  {tt)i , . . . ,  Wr^i }  be  a  set  of  arbitrary  positive  weights.  By  Taylor  expansion,  we  get 

r+l  r+1  r+1 


(3.5)  ^ «)«/(£>  MO  =  ^Wi/(£),M0  AI>\  +0(i''). 

^—1  i—l  i—l  \  / 

It  follows  from  (3.5)  and  (3.3)  that 

^.+1  ^  r+l  ^  /  df  ^  \  . 

(3.6)  mf{D  +  tAD,  Mj)  =  '^Wif{D,Mi)  -j-twr^i  /  —  (D,  Af^+i),  AD  \  -j-0{t^), 

i=i  i=i  \  ^  / 

By  (3.4),  we  can  choose  t  such  that  |^|  is  sufficiently  small  and 

<0. 
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Then  (3.6)  implies 


r+l  r+1 

Y,Wif{D  +  tAD,Mi}  <  Y,mf{D,Mi}. 

i=l  i=l 

That  is,  after  adding  one  more  design  point,  D  is  not  an  optimal  solution  no  matter  how  the  weights  are 
selected.  ■ 

In  summary,  if  the  number  of  design  points  is  less  than  the  number  of  free-design  variables,  then  it  is 
possible  to  have  a  first-order  reduction  of  the  value  of  /  at  some  Mach  number  Mj.^i  (in  order  of  |t|)  and  only 
second-order  increments  of  the  values  of  /  at  Af/s  (1  <i  <r)  (in  order  of  This  results  in  a  much  better 
performance  of  the  airfoil  at  an  off-design  point  with  a  marginal  deterioration  of  the  performance  at 

other  Mach  numbers. 

4,  Robust  Airfoil  Optimization  Formulations.  For  engineering  design  problems,  an  optimal  design 
is  usually  obtained  under  some  explicit /implicit  assumptions.  This  leads  to  a  design  that  works  well  under 
ideal  operating  conditions  but  may  perform  inadequately  under  non-ideal  (i.e.,  off- design)  conditions.  The 
problem  is  that  the  optimal  design  does  not  consider  the  uncertainty  or  variability  of  some  parameters/data 
that  will  affect  the  actual  performance  of  the  design  in  a  real-world  situation.  Therefore,  it  is  necessary  to 
include  uncertainties  in  a  practical  design  optimization  process.  In  this  paper,  we  assume  that  M  is  the  only 
uncertain  parameter  and  [Afinin?  ATmax]  is  the  range  of  Mach  numbers  considered. 

Let  p{M)  be  the  probability  density  function  of  the  uncertain  parameter  Af .  Then  a  stochastic  pro¬ 
gramming  formulation  for  airfoil  optimization  under  uncertainty  can  be  described  as  follows  [12]: 

(4.1)  min  /  c^{D,M,a{M))p{M)  dM 

DM') 

subject  to 

(4.2)  ci{D,  Af,  a(Af ))  >  c/  for  Afniin  <M<  Afniax 

and  D  ^  where  T  is  the  feasible  geometric  design  space  and  is  the  minimal  lift  value.  This  corresponds 
to  the  third  definition  of  ‘‘robust  optimization^^  given  in  the  introduction. 

On  the  other  hand^  one  can  also  use  the  following  rninirnax  optimization  formulation  for  robust  opti¬ 
mization  under  uncertainty  [4]: 

(4.3)  mill  max  p{M)ca{D,M,a{M)) 

subject  to  the  constraints  given  in  (4.2).  Here  p{M)  >  0  is  a  positive  weighting  function  of  Af.  This 
corresponds  to  the  second  definition  of  “robust  optimization”  given  in  the  introduction. 

The  constraints  on  the  lift  can  be  eliminated  by  choosing  an  appropriate  value  for  the  angle  of  attack 
corresponding  to  each  Af.  In  fact,  Elliott  and  Peraire  [6]  suggested  to  incorporate  a  means  for  adjusting 
the  angle  of  attack  to  satisfy  the  lift  constraint  into  the  flow  analysis  algorithm.  One  problem  with  this 
elimination  approach  is  that  the  angle  of  attack  becomes  a  function  of  Af  and  I>,  and  the  derivatives  of  a 
with  respect  to  D  have  to  be  computed  in  order  to  get  the  derivatives  of  c;  and  with  respect  to  D,  In 
this  study,  we  retain  the  lift  constraint  to  make  our  methods  applicable  to  general  constrained  aerodynamic 
optimization  problems  (with  constraints  on  the  pitching  moment  and  the  leading  edge  radius,  etc.). 
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5.  Approximations  to  Robust  Optimization  Formulations.  Since  we  cannot  compute  c;  and  q 
for  all  M  in  the  Mach  range  [Afmin,  Aimax]?  computationally  tractable  approximations  of  (4.1)  and  (4.3)  must 
be  used.  The  simplest  approximation  scheme  is  to  replace  [Afmin,  Afinax]  by  a  finite  subset  of  [Afmin,  Afinax]? 
say  {Ml,  M2, . . . ,  Mj.},  Then  (4.1)  is  reduced  to  the  following  multipoint  optimization  problem  [5]: 

r 

(5.1)  min  V  Wi  q(Z>,  Mi,  a*) 

i=l 

subject  to 

(5.2)  D  e  ci{D,Mi,ai)  >  cl  for  1  <  i  <  r, 

where  the  weights  wi  depend  on  the  probability  density  function  p{M)  and  a  chosen  integration  scheme. 
Similarly,  the  rninirnax  formulation  (4.3)  can  be  discretized  as  follows: 


(5.3) 


mill  max  pi  C(i{D,  Mi ,  ) , 


subject  to  the  constraints  given  in  (5.2).  Here  >  0  is  determined  by  p{M)  and  Mi. 

It  seems  that  (5.1)  and  (5.3)  are  completely  different.  However,  under  the  strict  complementarity 
condition  (i.e.,  Lagrange  multipliers  are  nonzero  for  all  active  constraints),  (5.1)  is  mathematically  equivalent 
to  (5.3).  Here  we  give  a  proof  of  the  mathematical  equivalence  between  (5.1)  and  (5.3). 

Let  (£>,  di, . . . ,  d^.)  be  a  stationary  solution  of  (5.1).  For  simplicity,  assume  that  D  is  in  the  interior  of  T 
and  the  equality  holds  for  lift  constraints.  Then,  under  the  strict  complementarity  condition,  the  following 
first-order  optimality  condition  (or  the  KKT-condition)  for  (5.1)  holds: 


(5.4) 


^Wi^{D,Mi,ai)  -  ^Xi^0,Mi,ai)  ^0, 
i—1  i—1 

Wi^^{b,Mi,ai)  —  \i^^{b,Mi,ai)  for  1  <  *  <  r, 
aai  aai 

Ai  >  0,  ci{D,  Mi,  d^)  —  cl  for  1  <  i  <  r. 


Define 


-y  —  ^  WiCdib,  Mi,  ai), 

i^l 


Pi  := 


7 

Cd{D,Mi,ai) 


>0 


for  1  <  i  <  r. 


6^^  :=  —  for  1  <  i  <  r. 
pi 


Then  (5.4)  implies  the  following  conditions: 


j2^i^{b,Mi,ai)^d, 

i—1  ^ 


(5.5) 


^ipi-^{£>,^^i,^i)  =  Xi-^{D,Mi,ai)  for  1  <  i  <  r, 
aai  aai 

Xi  >  0,  ci{D,  Mi,  di)  =  cl  for  1  <  i  <  r. 
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piCd{D,  Mi,  cci)  =j  for  1  <  i  <  r, 


9i  >  0  for  1  <  i  <  r,  and  ^  9i  =  1. 

i=l 

By  (5.5),  {j,D,di, , , , ,  dr)  satisfies  the  KKT-conditiori  of  the  following  optimization  problem: 

ry 

(5.6)  min  7  subject  to  ci{D,ai,  Mi)  >  Ci  and  Cd{D,ai,  Mi)  <  —  for  1  <  i  <  r. 

D,ai,j  pi 

Since  (5.6)  is  mathematically  equivalent  to  (5.3),  {D,  di, . . .  ,0:^.)  is  also  a  stationary  point  of  (5.3). 

On  the  other  hand,  assume  that  (£>,  di, . . . ,  dr)  is  a  stationary  point  of  (5.3),  the  equality  holds  in  (5.2), 
and  piCd{D,  Mi,di)  =  7  for  1  <i  <  r,  where 

7  max  Cd{D^Mi,di), 

Then  (7,  di, . . . ,  dr)  is  a  stationary  solution  of  (5.6).  If  we  further  assume  that  the  strict  complementarity 
condition  holds,  then  the  KKT-condition  (5.5)  for  (5.6)  holds.  It  follows  from  (5.5)  that  (5.4)  holds  for 
Wi  =  pi9i.  Therefore,  (£>, di, . . . ,  dr)  is  also  a  stationary  point  of  (5.1).  ■ 

Remark  1.  Note  that  the  strict  complementarity  condition  holds  naturally  for  (5.1).  In  fact,  if  the 
lift  is  greater  than  its  target  value,  then  one  can  reduce  the  angle  of  attack  so  that  the  lift  with  respect  to 
the  new  angle  of  attack  is  equal  to  the  target  value,  but  the  drag  with  respect  to  the  new  angle  of  attack 
is  reduced.  This  is  not  possible  at  a  stationary  point.  Therefore,  all  lift  constraints  at  a  stationary  point 
become  equality  constraints.  Also,  the  derivative  of  the  drag  with  respect  to  the  angle  of  attack  is  positive  for 
the  range  of  the  angle  of  attack  considered  herein,  which  implies  that  the  Lagrange  multipliers  are  positive 
(cf.  (5.4)  or  (5.5)).  Thus,  the  strict  complementarity  condition  for  (5.1)  always  holds.  An  implication  is 
that  it  is  possible  to  recover  an  optimal  solution  of  (5.1)  by  solving  (5.3)  with  appropriate  choices  of  pi^s. 
However,  it  is  a  little  more  difficult  to  claim  that  a  stationary  point  of  (5.3)  is  also  a  stationary  point 
of  (5.1).  In  the  proof  given  above,  we  used  two  additional  assumptions:  (a)  all  piCd{D,  AIi,di)  have  the 
same  value  7,  and  (b)  the  strict  complementarity  condition  holds  for  the  drag  constraints  in  (5.6).  The  first 
assumption  is  not  realistic  in  the  sense  that  we  do  not  know  how  to  choose  pi^s  to  make  piCd{D^  Mi,di)  have 
the  same  value  for  an  optimal  solution.  Adaptive  adjustments  of  piS  are  used  in  the  profile  optimization 
method  to  ensure  all  piCd{D,  Mi,  a^)  (1  <  i  <  r)  have  the  same  value  during  the  optimization  iterations.  The 
assumption  (b)  is  a  natural  one,  since  the  strict  complementarity  condition  for  the  drag  constraints  means 
that  every  drag  constraint  is  necessary  for  getting  the  optimal  solution.  ■ 

The  results  in  [5]  indicate  that  the  multipoint  airfoil  optimization  formulation  (5.1)  tends  to  produce 
airfoils  that  have  possibly  severe  degradation  in  the  off-design  performance,  which  Drela  referred  to  as  point- 
optimization  behavior.  The  above  equivalence  analysis  also  implies  that  the  same  conclusion  holds  for  the 
rninirnax  formulation  (5.3). 

Huyse  and  Lewis  [12]  concluded  that  the  point-optimization  behavior  can  be  attributed  to  the  discretiza¬ 
tion  error  between  (4.1)  and  (5.1).  The  multipoint  formulation  (5.1)  approximates  the  integral  as  a  discrete 
sum:  low  drag  is  obtained  at  the  selected  Mach  numbers  Afi, . . . ,  Af^.,  but  there  is  no  control  over,  nor 
requirements  for,  the  drag  at  the  other  Mach  numbers.  Due  to  the  highly  nonlinear  nature  of  the  PDEs 
for  flow  simulation,  the  optimizer  is  able  to  mold  the  objective  function  Cd  to  its  own  advantage:  distinct 
drag  troughs  appear  at  each  of  the  Mach  numbers  Mi, . . . ,  Mr  (see  also  [5]).  To  prevent  the  optimization 
algorithm  from  exploiting  the  approximation  error,  Huyse  and  Lewis  [12]  suggest  changing  the  integration 
points  Ml,...,  Mr  during  each  optimization  step.  This  effectively  ensures  that  low  drag  is  obtained  not 
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for  just  r  fixed  Mach  numbers,  but  for  any  combination  of  r  Mach  numbers  Afi, . . . ,  In  their  study, 
they  added  a  random  perturbation  to  the  integration  points  in  each  optimization  step  but  state  that  ‘‘any 
adaptive  scheme  that  varies  the  location  of  the  integration  points  Mk  for  each  optimization  step  will  do.” 

We  will  show  that  with  only  far  fewer  design  points  than  the  number  of  free-design  variables,  the  profile 
optimization  method  can  adaptively  adjust  p^’s  in  (5.3)  to  achieve  a  consistent  drag  reduction  over  the  given 
Mach  range,  so  there  will  be  no  degradation  in  the  off-design  performance  of  an  optimal  airfoil. 

6.  Profile  Optimization  Method,  In  this  section,  we  first  give  a  motivation  for  a  new  shape  opti¬ 
mization  strategy  based  on  the  robust  optimization  model  (5.3)  and  then  present  the  profile  optimization 
method. 

Let  D  be  a  fixed  feasible  design  vector.  We  can  plot  the  drag  with  respect  to  the  Mach  numbers  over 
the  interval  [Aimiiu  Afniax]  while  keeping  the  lift  at  a  constant  value  by  adjusting  the  angles  of  attack  for 
different  Mach  numbers.  Such  a  plot  will  be  called  a  profile  for  D  (see  Fig.  8.3(a)  for  example).  Ideally,  we 
want  to  obtain  a  design  vector  D  with  a  desirable  profile  (for  example,  with  the  drag  rise  occurring  at  a  very 
high  Mach  number) . 

Suppose  that  a  “perfect”  airfoil  D  could  be  obtained  by  some  magic  process,  which  might  be  either 
theoretical  or  experimental.  Let  us  look  at  the  unknown  process  from  a  reverse  engineering  point  of  view. 
That  is,  even  though  we  have  no  idea  about  how  the  airfoil  was  designed,  we  want  to  construct  an  optimization 
scheme  that  produces  an  airfoil  with  a  similar  or  better  profile. 

Assume  that  the  profile  curve  is  generated  by  computing  the  drag  coefficients  at  Mach  numbers  Afi,  Af2j 
. . . ,  Mr,  Let  1/ pi  be  the  value  of  the  drag  coefficient  of  the  “perfect”  airfoil  at  Mach  number  Ali  with  angle 
of  attack  a^,  which  is  chosen  so  that  the  lift  coefficient  is  equal  to  Then  the  global  optimal  value  of  (5.3) 
is  no  more  than  1.  If  (l),  di, . . . ,  d^.)  is  a  global  optimal  solution  of  (5.3),  then 

CdiDy  Mi,  ai)  <  —  =  Cd{D,  Mi,  ai)  for  1  <i  <r. 

Pi 

Therefore,  with  the  given  analysis  (ie.,  using  the  profile  of  the  drag  coefficient  at  Afi, . . . ,  Af^-),  D  is  no  worse 
than  the  “perfect”  design  vector  D. 

This  shows  that  for  any  given  profile  analysis  method,  it  is  possible  to  recover  or  improve  a  desirable 
airfoil  by  solving  (5.3),  if  the  weights  are  appropriately  chosen.  In  this  sense,  (5.3)  provides  a  very  flexible 
tool  for  airfoil  optimization.  By  experimenting  with  various  choices  of  the  weights  pi,  a  robust  airfoil  may 
be  obtained  (if  such  an  airfoil  exists).  However,  instead  of  guessing  the  appropriate  weights  in  (5.3),  we 
can  adaptively  adjust  the  weights  during  the  optimization  iterations.  This  leads  to  the  profile  optimization 
method. 

Before  presenting  the  profile  optimization  method,  we  give  a  linear  programming  formulation  of  a  trust- 
region  method  for  solving  (5.3).  Using  linearization  of  the  nonlinear  functions  in  the  mathematically  equiv¬ 
alent  formulation  (5.6)  at  the  current  iteration  point  {D^ ,  . . . ,  a^-,^),  we  get  the  following  trust-region 

subproblem  for  (5.6): 


(6.1) 


mill  7  subject  to 

A£>,Aa.j;,7 


—(jid  <  ADi  <  did  for  1  <  i  <  m, 

-ai^k  <  <  anmx  -  oci^k  for  1  <  i  <  r. 


C(  (I>*  +  (I>*  Ml),  AD 


+  ^iD^ai,k,Mi)Aai>ci 


for  1  <  i  <  r. 
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/  dcd 

\dD 


{D\ai,k 


,AfO,AZ>y 


da  pi 


for  1  <  i  <  r, 


where  A£>  and  Aa^  are  the  increments  for  D  and  respectively,  (Ji  and  6  are  positive  constants  that  define 
the  trust  region  for  AI>,  and  amax  is  the  maximum  angle  of  attack  allowed. 

Algorithm  6.1.  (Profile  Optimization  Method)  Suppose  that  is  a  given  design  vector  and 
Afi,  Af2, . . . ,  Mr  are  the  design  points.  Let  0  <  rj  <  1  be  the  predicted  percentage  reduction  rate  of  the  drag  for 
the  trust-region  method  and  let  0  <  e  <  1  be  a  parameter  for  termination  of  the  algorithm.  Then  construct 
a  new  design  vector  as  follows, 

(6.1.0)  Initialize  the  angles  of  attack:  Find  c^2,o?  •  •  •  ?  c^r,o  ^'^oh  that  ci{D^ Mi)  =  for  1  < 
i  <  r;  and  let  k  —  t), 

(6.1.1)  Adjust  weights:  Let 


1 

■  Cd{D^,ai^k,Mi) 


for  1  <  i  <  r. 


(6.1.2)  Check  early  termination:  If  the  zero  vector  is  an  optimal  solution  of  (6,1)^  then  output  as 
an  optimal  solution  and  terminate  the  algorithm, 

(6.1.3)  Find  a  trust  region  for  the  predicted  percentage  reduction  of  the  drag:  Find  6  >  0  such 
that  the  optimal  objective  function  value  of  (6,1)  is  (1  —  r/). 

(6.1.4)  Solve  the  trust-region  subproblem:  Find  the  least  norm  solution  (AI>^,  Aai^^, . . . ,  Aa^.^^)  of 

(6.1). 

(6.1.5)  Generate  the  new  iterate:  Let  +  Aa^^^  for  1  <  i  <  r%  and  +  AD^, 

(6.1.6)  Check  heuristic  termination  conditions:  If 


(6.2)  max  piC(i{D^^^ ,ai  k+i,  Mi)  >  1  and 

r 

(6.3)  tk  <  C(i{D^,ai^f^,Mi), 

i^i 

where  is  the  accumulative  reduction  of  the  drag  at  the  design  points  defined  by 

r 

tk  (cd{D\ai,k,Mi)  -  Cd{D>^+\ai^k+uMi)), 

i^l 

then  output  as  an  optimal  solution  and  terminate  the  algorithm, 

(6.1.7)  Start  a  new  iteration:  Update  k  k  -\-l  and  go  back  to  (6.1.1). 

Remark  2.  The  adaptive  adjustment  of  weights  makes  it  possible  to  achieve  a  consistent  drag  reduction 
over  the  Mach  range  [A/inin?  A/max]*  Note  that  moving  along  the  descent  direction  AD^  gives  a  simultaneous 
drag  reduction  at  all  design  points  Afi, . . . ,  Mr  (at  least  if  a  small  step  is  used).  If  the  selected  design  points 
are  a  “fair  representation”  of  all  Mach  numbers  in  [Afmin,  Afmax]?  then  moving  along  AI>^  will  not  induce  any 
hidden  rise  of  the  drag  at  some  unselected  Mach  number.  This  leads  to  a  robust  optimization  method  that 
achieves  a  consistent  drag  reduction  over  the  given  Mach  range  from  iteration  to  iteration  (see  the  profiles 
of  the  drag  coefficient  for  iterates  generated  by  the  profile  optimization  method  in  Fig.  8.4(a)). 

In  contrast,  without  adaptive  weight  adjustments,  it  is  necessary  to  use  at  least  {rri  +  I)  design  points, 
where  rri  is  the  number  of  free-design  variables,  in  order  to  avoid  point-optimization  at  selected  design  points 
(cf.  Sections  3  and  5). 

Besides  the  adaptive  adjustment  of  weights,  there  are  three  features  of  the  profile  optimization  method 
that  are  not  standard  for  a  trust-region  method. 
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(a)  First,  the  size  of  the  trust  region  is  modified  to  achieve  a  given  predicted  percentage  reduction  rate  r/ 
of  the  drag.  This  approach  is  employed  to  avoid  the  following  two  problems,  which  are  incurred  in  practical 
applications  of  a  standard  trust-region  method:  (i)  an  appropriate  size  of  the  trust  region  is  unknown  without 
prior  experiences,  and  (ii)  if  each  iteration  takes  a  long  time  to  complete,  then  rejecting  a  new  iterate  might 
be  hard  to  accept  by  a  designer. 

For  airfoil  shape  optimization,  if  the  percentage  reduction  of  q  is  too  small,  then  the  reduction  could  be 
a  consequence  of  numerical  errors  and  is  less  reliable.  If  the  predicted  percentage  reduction  of  is  too  large, 
then  the  size  of  the  trust  region  is  also  relatively  large  and  the  predicted  reduction  can  not  be  trusted.  By 
choosing  an  appropriate  percentage  reduction  rate  for  each  iteration,  these  two  problems  might  be  avoided. 

Note  that  the  choice  of  the  predicted  percentage  reduction  rate  depends  on  a  user’s  knowledge  of  the 
accuracy  of  the  simulation  analysis  tool  involved,  the  time  allowed  to  generate  a  new  design,  and  the  expected 
overall  performance  improvement.  For  example,  if  the  relative  numerical  error  for  computing  is  0.5%, 
each  iteration  takes  3  hours,  a  designer  has  1  day  to  generate  a  new  design,  and  the  expected  improvement 
is  8%,  then  the  designer  could  set  rf  —  1%  so  that  each  iteration  will  produce  some  meaningful  improvement 
of  the  design  and  the  final  design  might  have  about  8%  improvement  over  the  original  design  after  1  day  (or 
8  iterations). 

(b)  The  second  non-standard  feature  is  that  we  use  the  least  norm  solution  of  the  linear  programming 
problem  (6.1),  which  can  be  obtained  by  solving  a  quadratic  perturbation  of  (6.1)  [14].  Our  implicit  assump¬ 
tion  here  is  that  the  original  airfoil  is  reasonable.  Therefore,  we  do  not  want  a  new  airfoil  to  deviate  too  much 
from  the  original  one  if  unnecessary.  By  using  the  least  norm  solution  of  (6.1),  we  intend  to  select  a  new 
airfoil  closest  to  the  original  one  while  achieving  a  predetermined  amount  of  improvement  in  performance. 

(c)  The  third  non-standard  feature  of  the  profile  optimization  method  is  the  termination  criterion,  which 

is  based  on  design  heuristics.  The  basic  idea  is  that  if  a  good  descent  direction  for  drag  reduction  at  all 
design  points  does  not  work  (ie.,  (6.2)  holds)  and  the  overall  percentage  reduction  is  insignificant  (ie.,  (6.3) 
holds),  then  there  is  no  need  to  continue  the  optimization  process.  ■ 

7.  Implementation  of  Profile  Optimization  Method.  The  main  implementation  issues  related  to 
the  profile  optimization  method  are  the  following: 

•  parallel  processing  of  the  function  and  gradient  evaluations  for  r  Mach  numbers, 

•  finding  initial  values  of  the  angles  of  attack  for  r  Mach  numbers, 

•  finding  the  size  r/  for  the  trust  region  in  an  iteration, 

•  computing  the  least  norm  solution  of  (6.1). 

In  the  next  four  subsections,  we  discuss  these  four  implementation  issues. 

7.1.  Parallel  Processing  of  Function  and  Gradient  Evaluations.  Approximate  values  of  c;  (or 
Cd)  can  be  obtained  with  one  flow  analysis  (i.e.,  solving  the  PDE  once)  and  the  gradient  of  c;  (or  q)  can 
be  obtained  with  one  sensitivity  analysis  (i.e.,  solving  one  adjoint  equation).  The  function  values  and  the 
gradients  of  q  (or  q)  at  different  Mach  numbers  can  be  calculated  independently.  Therefore,  if  2r  processors 
are  available  in  a  parallel  computing  environment,  the  walltime  for  function  and  gradient  evaluations  in  each 
iteration  of  the  profile  optimization  is  one  flow  analysis  and  one  sensitivity  analysis. 

It  is  desirable  to  have  a  parallel  wrapper  code  for  the  profile  optimization  method  to  run  several  instances 
of  a  CFD  code  without  requiring  any  modifications  to  the  code.  For  this  purpose,  we  store  multiple  copies 
of  an  executable  CFD  code  (FUN2D  [15]  for  the  current  implementation)  under  different  directories  (called 
hosting  directories)  in  a  file  server.  Each  copy  of  the  CFD  code  communicates  with  the  profile  optimization 
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method  by  writing  its  output  into  its  own  hosting  directory  and  the  profile  optimization  method  collects 
all  function  values  and  gradients  from  the  hosting  directories.  The  architecture  of  a  parallel  computing 
environment  (such  as  shared  memory  or  distributed  memory)  has  no  impact  on  the  implementation  since 
there  is  no  communication  between  any  two  different  processors. 

The  current  parallel  wrapper  code  is  based  on  a  shell-level  parallel  computing  protocol  ‘‘pbsdsh”  on  Coral 
[13]  at  ICASE.  Coral  is  a  Beowulf-class  cluster  of  96  Intel  Pentium  processors  with  3  additional  file  servers. 
To  get  the  function  values  and  the  gradients  for  c;  and  at  r  different  Mach  numbers  simultaneously,  we 
just  create  2r  hosting  directories.  For  1  <  i  <  r,  in  the  i-th  directory,  we  use  FUN2D  to  compute  the 
function  value  and  the  gradient  of  c;  for  the  i-th  Mach  number  Af^,  and  in  the  (r  +  i)-th  directory,  we  use 
FUN  2D  to  compute  the  function  value  and  the  gradient  of  for  the  i-th  Mach  number  Mi,  Therefore,  with 
2r  processors,  we  get  all  function  values  and  gradients  for  ci  and  in  the  walltime  for  one  flow  analysis  and 
one  sensitivity  analysis. 

Note  that  on  Coral,  the  Unix  command  “pbsdsh  gradfun”  makes  the  same  executable  code  “gradfun” 
run  on  different  nodes  of  Coral  in  parallel.  Here  is  a  brief  description  of  how  to  write  a  C-code  that  generates 
an  executable  code  ‘^radfun”  for  parallel  processing  of  function  and  gradient  evaluations  on  Coral; 

•  Create  a  shell  script  called  “getdd” ,  whose  content  is 

echo  $PBD_NODENUM 

The  output  of  ‘‘getJd”  is  the  value  of  the  environment  variable  $PBDJ>^ODENUM,  which  is  (i  —  1) 
if  “get_id”  is  executed  on  the  i-th  node. 

•  Write  a  C-code  that  does  the  following: 

—  use  “gethostname( HOST-NAME)”  to  get  the  host  name  (such  as  n013)  of  the  node  running 
the  C-code; 

—  execute  a  system  command 

“get id  >  HOST-NAME” 

in  the  C-code  that  gets  the  $PBD-NODENUM  environment  variable  and  stores  it  in  a  file 
whose  name  is  the  content  of  the  string  HO  ST -NAME; 

—  read  the  index  number  id  in  the  file  whose  name  is  the  content  of  HO  ST -NAME; 

—  go  to  the  id-th  hosting  directory;  and 

—  execute  FUN2D  to  get  the  function  value  and  gradient  for  c;  (or  c^)  at  Mi  if  and  only  if  i  —id-\-l 
(or  r  +  i  =  id  +  1)  for  1  <  i  <  r. 


7,2,  Finding  Initial  Values  of  the  Angles  of  Attack.  For  a  fixed  airfoil  shape  and  a  fixed  Mach 
number,  the  lift  is  almost  a  linear  function  of  the  angle  of  attack  [10,  p.  99].  Therefore,  we  use  the  following 
search  method  based  on  a  linear  interpolation  to  find  an  angle  of  attack  for  which  the  corresponding  lift  is 
the  given  target  value. 

Algorithm  7.1.  Choose  an  error  tolerance  e  >  0  and  an  initial  adjustment  rate  0  <  k  <  1,  Let  ao  be 
an  initial  guess  of  the  angle  of  attack  for  a  given  Mach  number  M  and  a  given  geometric  design  vector  D. 
Then  the  following  procedure  will  find  a  value  of  a  such  that  |c;(I>,  Af,  a)  —  <  e. 

(7.1.0)  Initialize  the  range  variables  ai  and  a2:  If  ci{D^  Al^ao)  >  ^  then  ai  =  (1  — «;)ao  and  a2  =  ccq; 

otherwise^  ai  =  ao  and  a2  =  (1  +  Av)ao. 

(7.1.1)  Use  linear  interpolation  to  predict  a:  Let 


a  =  ai  + 


-  q(I>,  Af,ai) 
c;(I>,  Af,  a2)  -  c;(I>,  Af,  ai) 


(a2  -ai). 
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(7.1.2)  Check  the  termination  condition;  If 


\ci{D,M,a)  -c^l  <  ecf, 


then  output  a  and  terminate  the  algorithm, 

(7.1.3)  Update  ai  and  a2: 

•  if  a  <  ai^  then  a2  =  cci  and  ai  =  a; 

•  if  a  >  a2^  then  ai  =  a2  and  a2  =  cc; 

•  if  cci  <  a  <  a2  and  ci{D^  a)  <  then  ai  =  a; 

•  if  ai  <  a  <  a2  and  ci{D^  a)  >  then  a2  =  cc, 

(7.1.4)  Start  a  new  iteration:  Go  back  to  (7.1.1), 


7.3.  Finding  the  Size  of  the  Trust  Region.  It  is  well-known  that  the  optimal  objective  function 
value  of  (7.1)  is  a  non-increasing  and  piecewise  affine  function  of  6,  Therefore j  we  use  the  following  search 
method  based  on  a  linear  interpolation  to  find  6, 

Algorithm  7.2.  Let  'j{d)  be  the  optimal  function  value  of  (7,1)  for  any  given  6,  Choose  an  error 
tolerance  e  >  0  and  an  initial  adjustment  rate  0  <  k  <  1,  If  the  zero  vector  is  not  an  optimal  solution  of 

(7.1)  for  6  —  So  with  0  <  <  1?  then  for  any  t  >  0^  the  following  procedure  will  find  0  <  ^  <  1  such  that 

|7((J)  -(l-r/)|  <e. 

(7.2.0)  Initialize  the  range  variables  and  ^2^  7f  j (So)  >  1  —  then  Si  =  Jq  (^'i^d  S2  =  (1  +  k)So?' 
otherwise^  Si  =  (l  —  f^^So  and  S2  = 

(7.2.1)  Use  linear  interpolation  to  predict  S:  Let 


S  =  + 


(l-r/)  -7((5i) 

lih)  -  lih) 


Si). 


(7.2.2)  Check  the  termination  condition;  If 


ItW  -(!-»/)  I  <  e, 

then  output  S  and  terminate  the  algorithm, 

(7.2.3)  Update  Si  and  (!>2: 

•  if  S  <  Si  ^  then  S2  —  Si  and  Si  —  S; 

•  if  S  >  S2p  then  Si  =  S2  and  S2  =  S; 

•  if  Si  <  S  <  S2  and  j{S)  >  (1  —  r/)^  then  Si  —  S; 

•  if  Si  <  S  <  S2  and  j{S)  <  (1  —  r/)^  then  S2  =  S, 

(7.2.4)  Start  a  new  iteration:  Go  back  to  (7.2.1), 

Remark  3.  In  general,  we  use  the  value  of  S  in  the  previous  iteration  of  the  profile  optimization 
method  as  Jo*  With  this  choice  of  J,  ^  =  0.02,  and  e  =  0.01  •  r/,  we  usually  solve  about  4  linear  programming 
subproblems  to  find  S  in  our  numerical  simulation  runs.  ■ 

7.4.  Least  Norm  Solution  of  LP  Subproblem.  By  using  slack  variables  Ui  and  we  can  convert 
(6.1)  into  a  linear  programming  problem  with  equality  constraints  and  simple  bound  constraints: 


(7.1)  min  7  subject  to 

AD,Aai,j,Ui,Vi 

0  <Ui  <  C;*,  0  <  Ut  <  1,  for  1  <  i  <  r, 

<  aniax  -  for  1  <  i  <  r. 
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(jid  <  AA  <  did  for  1  <  i  <  : 


ct{D^,ai,k,Mi}  +  =  cj  for  1  <  i  <  r, 

c<^(I>^a^.*,A^^)  +  (^^{D^,ai,k,Mi},AD^  +  ^{D^,ai,k,Mi}Aai +Vi  =  J  for  1  <  i  <  r. 

The  least  norm  solution  of  (7.1)  can  be  found  by  solving  a  strictly  convex  quadratic  programming 
problem  [14]  that  is  derived  by  replacing  the  objective  function  in  (7.1)  by 


(7.2) 


J  +  -  I  7^  +  ^  [uf  +  u?  +  (Aa^)^]  +  ^2 


where  s  (=  10“^  in  our  implementation)  is  a  small  positive  scalar. 

Note  that  minimizing  the  objective  function  in  (7.2)  will  force  all  u^’s  to  be  more  or  less  the  same. 
As  a  consequence,  it  produces  a  descent  direction  that  gives  an  almost  identical  relative  reduction  of  the 
drag  for  all  the  design  points  Afi, . . . ,  Af^..  We  believe  that  this  balanced  reduction  of  drag  will  prevent 
point-optimization  at  the  sampled  design  points.  Furthermore,  the  least  norm  solution  generates  a  new 
design  vector  that  stays  as  close  to  the  original  design  vector  as  possible  while  achieving  the  predetermined 
percentage  reduction  of  the  drag.  This  is  desirable  in  a  practical  design  process,  since  if  a  new  design  deviates 
too  far  away  from  the  original  airfoil,  then  designers  might  have  less  confidence  in  the  predicted  performance 
of  the  new  airfoil. 


8.  Numerical  Simulation  Results.  We  use  the  NACA-0012  airfoil  as  the  initial  point  for  all  opti¬ 
mization  methods  discussed  in  this  section.  A  periodic  spline  representation  of  the  NACA-0012  with  23 
control  points  is  used  to  get  an  initial  parametric  model  of  the  airfoil: 

22  22 

x^'^aiPi{t),  for  0  <  f  <  1, 

where  pi  and  qi  are  spline  basis  functions.  The  shape  of  the  airfoil  can  be  modified  by  changing  the  values 
of  the  The  tip  and  the  tail  of  the  airfoil  are  fixed  during  the  optimization  (ie.,  6o,  ^n,  and  ^>22 

are  fixed)  so  that  the  chord  length  remains  constant.  As  a  consequence,  we  have  20  free-design  variables 
(^1, . . . ,  ^12?  •  •  •  j  ^21)  iu  the  airfoil  shape  optimization. 

The  CFD  code  FUN2D  is  used  to  compute  the  function  values  and  gradients  of  the  lift  and  drag  [15]. 
Some  technical  details  of  the  flow  solver  (for  function  evaluations)  and  the  adjoint  equation  solver  (for 
gradient  evaluations)  in  FUN2D  can  be  found  in  [1]  and  [16],  respectively. 

We  use  Euler  flow  to  demonstrate  the  usefulness  of  the  profile  optimization  method  as  a  robust  airfoil 
optimization  tool.  The  main  reason  to  choose  Euler  flow  simulation  analysis  is  that  we  can  complete  our 
preliminary  evaluation  of  the  profile  optimization  method  in  a  reasonable  amount  of  time. 

Elliott  and  Peraire  [6,  Section  1]  stated  that  almost  any  drag  minimization  exercise  based  on  the  Euler 
equations  and  applied  to  modern  supercritical  wings  in  cruise  condition  is  doomed  to  failure.  However,  they 
also  concluded  that  Euler-based  optimization  can  be  useful  as  a  first  step  in  the  design  process  [7,  Abstract]. 

For  a  given  initial  grid,  FUN2D  generates  numerical  approximations  to  the  lift /drag  coefficients  and 
their  derivatives.  We  use  10“^  as  the  error  tolerance  for  the  flow  solver  in  FUN2D  to  find  an  approximate 
solution  of  the  nonlinear  system  derived  from  discretization  of  the  Euler  equation.  The  number  of  time  steps 
used  in  solving  the  discrete  adjoint  equation  is  10000  and  the  spring  analogy  method  in  FUN2D  is  used  for 
grid  movement  [16]. 
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The  initial  grid  for  FUN2D  used  in  this  paper  has  124  points  around  the  airfoil  and  32  points  at  the  far 
field  (at  a  distance  of  50  chord  lengths).  The  grid  has  3060  nodes,  9025  faces,  and  6120  elements.  Fig.  8.1 
shows  the  initial  grid  around  the  NACA-0012.  See  [11]  for  a  comparison  of  the  numerical  approximations  of 
the  lift  and  drag  using  this  grid  and  other  grids. 


Fig.  8,1.  The  grid  used  for  solving  the  Euler  equation  by  FUN2D 

All  simulation  results  are  obtained  by  parallel  processing  of  function  and  gradient  evaluations  on  Coral 
(cf.  Subsection  7.1).  For  the  grid  we  use,  it  takes  about  120  seconds  to  solve  a  flow  problem  and  about  180 
seconds  to  solve  an  adjoint  equation.  Therefore,  with  2r  processors,  each  iteration  of  the  profile  optimization 
method  takes  about  5  minutes,  no  matter  what  the  number  of  design  points  is. 

We  compute  the  least  norm  solution  of  (7.1)  by  the  quadratic  programming  solver  QLD,  which  was  first 
developed  by  Powell  [18]  and  then  modified  by  Schittkowski  [19]. 

Our  numerical  simulations  are  designed  to  examine  the  impact  of  the  number  of  design  points  and 
the  target  lift  value  on  the  profile  optimization  method.  Also,  we  shall  compare  the  simulation  results 
obtained  by  the  profile  optimization  method  with  the  multipoint  optimization  method  and  the  expected 
value  optimization  method.  Therefore,  we  consider  the  following  cases  for  numerical  simulation: 

1.  the  profile  optimization  with  4  design  points  equally  spaced  in  the  Mach  range  [0.7, 0.8],  =  0.4,  a 

fixed  percentage  reduction  rate  r/  =  5%,  and  a  termination  parameter  e  =  0.3%; 

2.  the  profile  optimization  with  3  design  points  equally  spaced  in  the  Mach  range  [0.7, 0.8],  =  0.4,  a 

fixed  percentage  reduction  rate  r/  =  5%,  and  a  termination  parameter  e  =  0.3%; 

3.  the  profile  optimization  with  8  design  points  equally  spaced  in  the  Mach  range  [0.7, 0.8],  =  0.4,  a 

fixed  percentage  reduction  rate  r/  =  5%,  and  a  termination  parameter  e  =  0.3%; 

4.  the  profile  optimization  with  4  design  points  equally  spaced  in  the  Mach  range  [0.7, 0.8],  =0.2,  a 

fixed  percentage  reduction  rate  r/  =  5%,  and  a  termination  parameter  e  =  0.3%; 

5.  the  multipoint  airfoil  optimization  with  4  design  points  equally  spaced  in  the  Mach  range  [0.7, 0.8], 

=  1/6,  =  1/3,  and  =  0.4; 

6.  the  expected  value  airfoil  optimization  with  adaptive  changes  of  Wi  corresponding  to  4  randomized 
integration  points,  and  =0.4. 

The  first  case  is  the  benchmark  case.  The  next  two  cases  are  included  to  examine  the  impact  of  the 
number  of  design  points  on  the  profile  optimization.  We  use  Case  4  to  examine  the  impact  of  the  target 
lift  value  on  the  profile  optimization.  The  last  two  cases  are  for  comparison  of  the  profile  optimization  with 
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other  optimization  methods.  The  profiles  of  the  drag  coefficient  for  a  given  are  plotted  by  using  24  equally 
spaced  Mach  numbers  in  [0.7, 0.8]. 

8.1.  Impact  of  the  Number  of  Design  Points.  The  profile  optimization  method  terminates  after 
69,  67,  and  57  iterations  for  3,  4,  and  8  design  points,  respectively.  The  number  of  design  points  has  no 
significant  impact  on  the  optimal  airfoils  generated  by  the  profile  optimization  method.  All  optimal  airfoils 
have  similar  shapes  (cf.  Table  8.1)  and  only  a  marginal  drag  reduction  is  achieved  by  adding  more  iterations 
(cf.  Table  8.2). 


Table  8.1 

The  relative  differences  (in  percentage)  of  the  design  vectors 


NACA-0012 

3  Points 

4  Points 

8  Points 

NACA-0012 

- 

24.5 

25.7 

23.8 

3  Points 

24.7 

- 

1.7 

1.2 

4  Points 

25.9 

1.8 

- 

2.6 

8  Points 

23.9 

1.2 

2.6 

- 

Each  number  in  Table  8.1  denotes  the  relative  difference  of  two  design  vectors:  ||D^.  —  Dc||/||D^.||,  where 
D^.  (or  Dc)  is  either  the  NACA-0012  or  the  optimal  design  vector  obtained  by  the  profile  optimization 
method  with  the  number  of  design  points  listed  in  the  corresponding  row  (or  column) . 

A  typical  optimal  airfoil  is  given  in  Fig.  8.2(b).  This  Mach  wave  plot  illustrates  how  one  strong  shock 
wave  for  the  NACA-0012  is  reduced  to  several  weaker  shock  waves  for  the  optimal  airfoil. 


(a)  NACA-00 1 2  (b)  Profile  Optimization 


Fig,  8.2.  The  Mach  waves  (for  the  Mach  number  0.796)  around  the  NACA-0012  and  an  optimal  airfoil  generated  by  the 
profile  optimization  method 

With  4  design  points,  the  profile  optimization  method  achieves  consistent  drag  reduction  over  the  Mach 
range  [0.7, 0.8]  (cf.  Fig.  8.3(a)  and  8.4(a)).  There  is  no  random  distortion  of  airfoil  shapes  during  the 
optimization  process  (cf.  Fig.  8.3(b)  and  8.4(b)). 

It  is  worth  pointing  out  that  the  iterates  generated  by  the  profile  optimization  method  are  also  indepen¬ 
dent  of  the  number  of  design  points.  For  example,  the  56-th  iterate  generated  by  the  profile  optimization 
method  for  4  design  points  differs  from  the  56-th  iterate  corresponding  to  3  (or  8)  design  points  by  0.35% 
(or  1.2%).  The  profiles  of  the  drag  coefficient  for  optimal  airfoils  generated  by  the  profile  optimization  with 
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3  and  4  design  points,  respectively,  are  almost  identical.  The  profile  of  the  drag  coefficient  for  the  optimal 
airfoil  corresponding  to  8  design  points  is  very  similar  to  the  drag  profile  of  the  57-th  iterate  generated  by 
the  profile  optimization  with  4  design  points. 

8.2.  Impact  of  the  Target  Lift  Value.  The  target  lift  value  has  a  quite  significant  impact  on  the 
optimal  airfoil  generated  by  the  profile  optimization  method. 

(b)  Profile  Optimization  With  Lift  at  0.2 

0.06 

0.04 

0.02 

0 

-0.02 

-0.04 

-0.06 

0.7  0.71  0.72  0.73  0.74  0.75  0.76  0.77  0.78  0.79  0.8 

Mach  Number  From  Left  to  Right;  NACA-0012,  6-th,  12-th,  18-th,  24-th,  and  33-th  Iterates 

Fig,  8,3.  The  two  figures  illustrate  the  behaviors  of  the  itemtes  generuted  by  the  profile  optimization  method^  with  cj  —  0,2 
and  4  equally  spaced  design  points.  Fig.  8.3(a)  shows  profiles  of  the  drag  coefficient  versus  the  Mach  number  for  the  itemtes 
and  Fig.  8.3(b)  shows  the  shapes  of  airfoils  corresponding  to  the  iterates.  (Note.  The  airfoils  are  shifted  to  improve  the 
visibility.) 

With  —  0.2,  the  NACA-0012  is  almost  optimal  for  low  Mach  numbers.  Therefore,  we  can  not 
simultaneously  reduce  the  drag  over  the  Mach  range  [0.7, 0.8].  The  profile  optimization  method  tries  to  keep 
the  drag  as  low  as  possible  for  Mach  numbers  between  0.7  and  0.75,  while  reducing  the  drag  significantly  for 
Mach  numbers  between  0.75  and  0.8  (cf.  Fig.  8.3(a)). 

Note  that  the  drag  bucket  (i.e.,  a  drop  of  the  drag  before  its  dramatic  rise)  occurs  near  M  —  0.73  for 
the  NACA-0012  and  the  drag  bucket  for  the  optimal  airfoil  occurs  near  M  =  0.78.  Such  a  delay  of  the  drag 
rise  is  desirable.  The  only  compromise  made  by  the  optimal  airfoil  to  the  NACA-0012  is  a  small  increase  of 
the  drag  around  the  original  drag  bucket  for  the  NACA-0012,  which  is  quite  reasonable. 

For  the  higher  target  lift  =  0.4,  the  optimal  airfoil  shape  (shown  in  Fig.  8.4(b))  deviates  more  from 
the  NACA-0012.  Also,  the  33rd  iterate  generated  by  the  profile  optimization  method  for  =0.4  differs 
from  the  33rd  iterate  for  =  0.2  by  10%. 

Note  that  the  shape  variations  of  the  airfoils  corresponding  to  iterates  generated  by  the  profile  optimiza¬ 
tion  method  follow  a  similar  trend  no  matter  what  the  target  lift  is:  the  aft  end  thickens  while  the  front 
section  narrows.  This  pushes  the  shock  location  towards  the  aft  region  of  the  airfoil. 

A  similar  airfoil  with  thin  front  and  fat  tail  was  obtained  by  Elliott  and  Peraire  [6,  Figure  34]  in  a 
totally  different  context,  where  they  used  two  thickness  design  variables  and  one  camber  design  variable  for 
a  lift-constrained  drag  optimization  with  2-D  separated  viscous  flow  and  an  additional  area  constraint.  The 
initial  airfoil  used  by  them  is  also  the  NACA-0012.  Their  explanation  of  the  merit  of  an  airfoil  with  thinner 
front  and  fatter  tail  over  the  NACA-0012  is  that  the  thickness  distribution  has  been  redistributed  such  that 
the  maximum  is  further  aft  (i.e.,  closer  to  the  tail),  like  the  early  natural  laminar  flow  airfoils,  this  delays 
the  start  of  the  adverse  pressure  gradient  to  aft  of  that  maximum  thickness  point. 
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8.3.  Impact  of  the  Optimization  Strategy.  From  the  analysis  given  in  Sections  3  and  5,  we  know 
that  if  the  design  points  and  weights  are  fixed  as  in  either  (5.1)  or  (5.3),  then  a  severe  degradation  in  the 
off-design  performance  is  likely  when  the  number  of  design  points  (i.e.,  4)  is  much  smaller  than  the  number 
of  free-design  variables  (i.e.,  20).  Fig.  8.5(a)  clearly  reveals  the  point-optimization  features  of  the  optimal 
airfoil  generated  by  the  multipoint  optimization  method.  Fig.  8.4(c)  suggests  that  the  point-optimization 
behavior  does  not  occur  until  the  end  of  the  optimization  iterations,  since  the  drag  curve  for  the  25-th  iterate 
does  not  have  any  drag  trough. 

The  expected  value  optimization  method  [12]  adds  a  random  perturbation  to  the  integration  points 
from  iteration  to  iteration  and  adjusts  the  weights  accordingly  so  that  a  severe  degradation  in  the  off-design 
performance  can  be  avoided.  Fig.  8.4(e)  shows  that  this  technique  allows  the  optimizer  to  reduce  the  drag 
while  avoiding  point-optimization. 

The  profile  optimization  avoids  the  problem  of  point-optimization  at  the  design  points  (cf.  Fig.  8.3(a) 
and  8.4(a))  by  using  a  smart  descent  direction  to  achieve  consistent  drag  reduction  over  the  whole  Mach 
range. 

Table  8,2 

Relative  drag  reduetion  rates  (in  percentage) 


Max 

Min 

Average 

Profile  Optimization  (3  Points) 

82.9 

5.1 

70.2 

Profile  Optimization  (4  Points) 

82.7 

1.4 

69.7 

Profile  Optimization  (8  Points) 

82.6 

4.1 

67.9 

Multipoint  Optimization  (4  Points) 

89.4 

46.2 

81.1 

Expected  Value  Optimization  (4  Points) 

89.4 

49.6 

82.6 

Table  8.2  gives  an  overall  sense  of  how  each  optimization  method  reduces  the  drag  coefficient  with 
c*  =  0.4.  Let  AfijAf2,...,  M24,  be  24  equally  spaced  points  in  [0.7, 0.8].  For  each  Mi^  let  (or  Cd^i)  be 
the  drag  coefficient  of  the  NACA-0012  (or  an  optimal  airfoil)  with  the  lift  coefficient  fixed  at  0.4.  Then  the 
relative  reduction  of  the  drag  at  each  Mach  number  AIi  is  {cd^i  —  Cd^i) ! The  numbers  in  “Max”^  “Min”, 
and  “Average”  columns  of  Table  8.2  correspond  to 


max 

l<i<24 


(^d,i  —  ^d,i 
^d,i 


rnin  — ^ ^ ,  and 

l<i<24  Cd^i 


ed,i  -  ^  j  Cd^  ,  respectively. 

Note  also  that  the  expected  value  optimization  method  (as  well  as  the  multipoint  optimization  method) 
resulted  in  an  optimal  airfoil  with  a  drag  curve  that  is  consistently  lower  over  the  Mach  range  than  the 
profile  optimization  method  (cf.  Fig.  8.5(a)).  It  seems  that  it  should  be  possible  for  the  profile  optimization 
method  to  find  a  descent  direction  for  a  consistent  drag  reduction  over  the  whole  Mach  range  [0.7, 0.8],  by 
moving  toward  the  optimal  solution  given  by  either  the  multipoint  optimization  method  or  the  expected 
value  optimization  method.  However,  due  to  nonlinearity  of  the  drag  coefficient  (with  respect  to  the  design 
variables),  the  profile  optimization  method  can  not  find  such  a  descent  direction  by  using  the  local  derivative 
information  only. 


9,  Conclusion,  In  this  paper,  we  introduce  a  new  robust  optimization  scheme  called  the  profile  opti¬ 
mization  method  and  use  airfoil  optimization  under  uncertain  flight  conditions  as  a  case  study  to  evaluate 
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(a)  Profile  Optimization 


(b)  Profile  Optimization 


From  Left  to  Right;  NACA-0012,  13-th.  26-th,  39-th.  52-th,  and  69-th  Iterates 


(c)  Multipoint  Optimization 


(d)  Multipoint  Optimization 


From  Left  to  Right;  NACA-0012,  6-th,  12-th,  18-th,  25-th.  and  33-th  Iterates 


(e)  Expected  Value  Optimization 


(f)  Expected  Value  Optimization 


Fig.  8.4.  These  figures  show  the  behaviors  of  iterates  generated  by  various  optimization  methods  with  cj  —  0.4  and  4 
design  points. 


this  method.  We  used  an  Euler-based  CFD  code,  which  does  not  include  viscous  effects,  to  test  the  profile 
optimization  method.  Because  of  this  lack  of  fidelity,  the  generated  airfoils  may  be  somewhat  unrealistic. 

The  profile  optimization  method  adaptively  adjusts  the  weights  in  a  rninirnax  optimization  formulation 
to  find  a  drag  reduction  direction  for  all  design  conditions,  which  leads  to  a  consistent  drag  reduction  over 
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^  ^q-3  (a)  Profiles  of  the  Drag  Coefficient 


(b)  Optimal  Airfoils 


From  Left  to  Right;  Profile  Opt.  Multipoint  Opt,  and  Expected  Value  Opt 


Fig.  8.5.  Fig.  8.5(a)  compares  the  profiles  of  the  drug  coefficient  for  various  optimal  airfoils  and  Fig.  8.5(b)  comjyares 
the  optimal  airfoil  shapes. 


a  given  range  of  Mach  numbers  from  iteration  to  iteration. 

Without  adaptive  adjustment  of  weights^  we  prove  that  it  is  necessary  to  use  at  least  (m  +  1)  design 
points j  where  rri  is  the  number  of  free-design  variables j  in  order  to  avoid  point-optimization  at  selected  design 
points. 

Our  numerical  results  demonstrate  that  the  profile  optimization  method  is  not  sensitive  to  the  number  of 
selected  design  points.  With  20  free  geometric  design  variables  and  as  little  as  4  design  points ^  the  optimal 
airfoil  generated  by  the  profile  optimization  method  is  smooth  and  has  no  degradation  in  the  off-design 
performance. 

The  profile  optimization  method  can  be  easily  modified  to  solve  other  optimization  problems  under 
uncertainty  by  replacing  q  with  another  performance  measurement  function  and  M  with  other  uncertain 
parameters. 

The  profile  optimization  method  also  has  the  potential  of  becoming  a  practical  design  tool  for  optimiza¬ 
tion  under  uncertainty.  The  use  of  a  small  number  of  sampled  design  points  from  the  range  of  uncertain 
parameters  makes  the  profile  optimization  method  computationally  affordable.  The  consistent  performance 
improvements  over  the  range  of  uncertain  parameters  from  iteration  to  iteration  allows  a  designer  to  stop 
the  iterative  process  at  any  time  with  an  improved  new  design. 
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